DNA Unzipping under a periodic force: Scaling behavior of hysteresis loop 
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We propose a dynamic transition, which has potential application in understanding non- 
equilibrium processes such as transcription and replication of DNA, unfolding of proteins and RNA 
etc. Here, we demonstrate that without changing the physiological condition, a DNA may be brought 
to a new dynamic state from the zipped or unzipped state by varying the frequency of the applied 
force. We report the value of various critical exponents associated with this transition, which are 
amenable to verification in the force spectroscopic experiments. 
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The mechanism involved in the separation of a double 
stranded DNA (dsDNA) into two single stranded DNA 
(ssDNA) is a prerequisite for understanding the pro- 
cesses like replication and transcription. In vitro, open- 
ing of DNA is generally achieved either by increasing the 
temperature (85-90°C) termed as thermal melting or by 
changing the pH value of the solvent (< 3 or > 9) called 
DNA denaturation |2|. However, such a drastic change 
in the physiological conditions is not possible in living 
systems. The mechanism of opening of dsDNA in vivo 
is quite complex and is initiated by helicases, DNA and 
RNA polymerase, ssb proteins etc, which exert force of 
the order of pN and as a result DNA unwinds. It is 
now possible to unzip the two strands of a DNA using 
techniques like optical tweezers, atomic force microscopy, 
magnetic tweezers etc [H, 0] • Theoretical understandings 
of DNA unzipping are mostly based on equilibrium con- 
dition 



However, living systems are open systems and never 
at equilibrium. Understanding the separation of DNA in 
the equilibrium is one approach, but another route is to 
perform the analysis in a situation, which closely resem- 
bles the living systems i. e. in non-equilibrium condition. 
Moreover, DNA helicases are ATP driven molecular mo- 
tors. The periodic hydrolysis of ATP to ADP can gen- 
erate continuously push and pull kind of motion, which 
spans a wide range of length and time scales. In the cat- 
alytic sites, it is localized (~ 1A) and has fast chemical 
reaction (~ fs), whereas, it can extend to the large scale 
(up to 10 nm) and has much slower conformational mo- 
tions (up to ms) in physiological functions [l|, [9j . This 
along with other studies suggest that the force acting on 
DNA is oscillatory in nature of varyin g fre quency (v) and 
amplitude (F) rather than constant |1CH12| |. The typi- 
cal time scale of fluctuations involved in unzipping varies 
from ns to (is, and therefore, it is neccessary to under- 
stand the influence of oscillatory force on the unzipping. 
Surprisingly, in most of the studies (experiments, theo- 
ries and simulations), the force or loading rate is kept 
constant (l3l |. 

In DNA unzipping, the equilibrium response of reac- 
tion coordinate (extension (y)) to the constant force is 
well understood in terms of simple models amenable to 
statistical mechanics (T3H15I . However, when a dsDNA 



is driven by an oscillatory force, y will also oscillate and 
lag behind the force due to the relaxation delay. This 
relaxation delay induces hysteresis in the force-extension 
(/ — y) curve of DNA unzipping, which has been recently 
observed in simulations and experiment [T6Hl8j . The na- 
ture of hysteresis and its dependence on the F and v is 
well studied in context of the spin systems [l9l - [2l| . It is 
found that the area under the hysteresis loop Ai oop scales 
as F a v^ . Values of cv and j3 differ from system to system 
[2ll ] . However, for DNA unzipping, the non-equilibrium 
response of y to the oscillatory force remains elusive. 

In this Letter, we shall probe how a living system not in 
equilibrium remains structurally stable for an extended 
period of time. More specifically, we show that under 
a certain physiological condition, a dsDNA remains in 
the steady and stable (zipped or open) state for an ex- 
tended period of time. Furthermore, without any change 
in temperature (T) or pH of the solvent, by varying v 
alone, a dsDNA may be brought from the time averaged 
open or the zipped state to a new dynamic state (hys- 
teretic), oscillating between the zipped and the unzipped 
states, which is dynamical in origin and vanishes in the 
quasi static limit [21| . We evaluate the scaling exponents 
a and j3 associated with Ai oop , which are amenable to 
verification in the force spectroscopic experiments. We 
also show that using the work theorem [221 ] . it is pos- 
sible to extract the equilibrium / — y curve from the 
non-equilibrium pathways. 

In order to study the dynamical stability of DNA un- 
der the oscillatory force, one may fix v and vary F or 
keep F constant and vary v. Since, time scale involved 
here is quite large, therefore, an all atom simulation of 
longer base sequence is beyond our computer resource. 
In view of this, we adopted minimal off-lattice coarse 
grained model of DNA. Here, a bead represents few bases 
associated with sugar and phosphate groups. We con- 
sider a sequence in such a way that the first half of the 
chain is complementary to the other half. This gives the 
possibility of the formation of a dsdNA at low tempera- 
ture [23| . The energy of the model system, we adopted, 
is defined as 
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FIG. 1: DNA in a zipped and open conformation. One end is 
fixed and the other end is subjected to a pulling force. 



where N(= 32) is 
between beads da. 



the number of beads, 
is defined as |r, — fA 



The distance 
where r, and 



fj denote the position of bead i and j, respectively. The 
harmonic (first) term with spring constant k (=100) cou- 
ples the adjacent beads along the chain. The remaining 
terms correspond to Lennard-Jones (LJ) potential. The 
first term of LJ potential takes care of the "excluded vol- 
ume effect", where we set B = 1. We assign the base 
pairing interaction Ag — 1 for native contacts and 
for non- native ones (l6| . This choice corresponds to the 
Go model [23|. By native, we mean that the first base 
forms pair with the N th (last one) base only and second 
base with (N — l) th base and so on as shown in Fig.l. 
The parameter do(= 1.12) corresponds to the equilibrium 
distance in the harmonic potential, which is close to the 
equilibrium position of the average L- J potential. In the 
Hamiltonian (Eq. 1), we use dimensionless distances and 
energy parameters. We obtained the dynamics by using 
the following Langevin equation [2£ 



dt 2 



^ dt 
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where m(= 1) and £(= 0.4) are the mass of a bead and 
friction coefficient, respectively. Here, F c is defined as 
— and the random force V is a white noise (2(| i.e., 
< f (t)r(i') >= 2(T5(t - t'). This keeps temperature 
constant throughout the simulation. The equation of 
motion is integrated using 6 th order predictor corrector 
algorithm with a time step St = 0.025 [26[. We add an 
energy —f(t).y(t) to the total energy of the system given 
by Eq.l The value of / increases to its maximum value 
F in the m s steps at the interval A/(= 0.01) and then 
decreases to in the same way. Since, we are interested 
in the non-equilibrium regime, we allow only n LD time 
steps (much below the equilibrium time) in each incre- 
ment of A/. Here, y(f) is the distance between two ends 
of the chain after the n LD time steps, where the force 
/ is applied. We keep sum of the time spent (t = 2nm s ) 
in each force cycle constant to keep v — (1/t) constant. 

We have recently shown that the equilibrium force- 
temperature diagram of a short DNA is in good agree- 
ment with the two state model in the entire range of the 
/ and T 11611 . In the following, we keep T = 0.1 and 
F > 0.25 27]. In Fig.2, we plot the value of y(f) (av- 
eraged over C =1000 cycles) with / for different values 
of v. It is interesting to note that the average value of 
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FIG. 2: (y) of DNA as a function of the cyclic force of ampli- 
tudes (a) 0.4 and (b) 1.0 at different v. (a) At high v, DNA 
remains in the zipped state with a small hysteresis loop. As 
v decreases, the system extends from the zipped state to the 
open state with a bigger loop. For v — > 0, the hysteresis loop 
vanishes and the system approaches to the equilibrium path, 
(b) DNA remains in the open state at high v and approaches 
the equilibrium path from below as v — > 0. 



y(f) for different initial conformations remains almost 
the same, showing that the system is in the steady state 
[28j ] . All plots show hysteresis, however, the dynamic or- 
der parameter Ai oop = § y.df depends upon F and v. 
The melting tempearture is usually defined, when half of 
the base-pairs of the chain gets open 0- If y(f) is less 
than 5, we call the system is in the zipped state, where 
as if y(f) > 5, it is in the open state or unzipped state. 
At high v, it is evident that for small F, dsDNA remains 
in the zipped state (Fig. 2a), whereas at high F, it is in 
the unzipped state (Fig. 2b), irrespective of initial con- 
formation. Moreover, the path of y(f) for the force to 
F is different than that of the path for the force F to 0, 
which constitutes a hysteresis loop. Decrease in v leads 
to a bigger path of the hysteresis loop. Depending on the 
amplitude, system starts from the zipped conformation 
as shown in Fig. 2a (or open conformations shown in Fig. 
2b) and then gradually approaches to the open state (or 
the zipped state) and back to the initial state. 

One may note that even though / decreases from F to 
(Fig. 2a), y(f) increases and there is some lag, after 
which it decreases. One may recall that the relaxation 
time is much higher compare to the time spent at each 
interval of A/, and hence, increase in y(f) with decrease 
in / indicates that system gets more time to relax. As 
a result y(f) approaches to a path, which is close to the 
equilibrium. Once the system gets enough time, the lag 
disappears. Similar lag one can expect, when the system 
starts from the open state at high v. However, in this 
case as v decreases, y(f) decreases with increase in / 
(Fig. 2b). In both cases , whether dsDNA starts from the 
zipped or open state, as v —> 0, the system approaches 
to the equilibrium / — x curve and Ai oop vanishes. 

The area under the loop is the measure of the energy 
dissipated over a cycle. Moreover, in both cases (Fig. 2a 
& 2b), Ai 00 p is quite small at high v. At fixed F, as v 
decreases, Ai oop increases (Fig. 3a). At some value of u, 
Aioop is found to be maximum and then it approaches to 
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FIG. 3: Variation of loop area Ai oop (a) with v at different 
F, (b) with F at different v. 



zero. In Fig. 3b, we have depicted the plots of Ai oop with 
F for different v. In this case also, Ai oop increases with 
F and above a certain amplitude, it starts decreasing. 
However, one can observe here that A is the slow variant 
of F and never approaches to zero for all F. 

The other dynamic order parameter Q = — § x(t)dt, 
studied in context of magnetic system, has recently been 
applied to obtain the the force-frequency (F — v) dia- 
gram of DNA hairpin (29|. The explicit loop of hairpin 
captures essential features of a long DNA with implicit 
bubbles, and / — T diagram is found to be qualitatively 
similar to the experiment [4J. In Fig. 4, we plot the 
value Q with cycles for different v and F. The distri- 
bution shows that the path remains in zipped state or 
in open state or in dynamic state, depending on F and 
v. In Contrast to DNA hairpin, a short DNA shows a 
continuous transition from the zipped state to the new 
dynamic state as frequency decreases. 

We now focus our study on the scaling of the area of 
hysteresis loop. In Fig. 5a, we have plotted Ai oop as a 
function of (i*V) - 5 . For low v, we observe that all plots 
for different F collapse on a straight line. This gives the 
value of a = 0.5 = /3, which are consistent with the mean 
field values for a time dependedent hysteretic response to 
periodic force in case of the isotropic spin 20]. At high v, 
depending on the amplitude, the system remains either 
in the zipped state (low F) or in the open state (high F) . 
In contrast to the spin system, in this case two states are 
asymmetric and hence, one should not expect the same 
scaling. In fact, at low F, we find that Ai oop scales as 
v^ 1 f 5 , which is different than the scaling predicted for 
the spin system, where it scales as v~ x f 2 [20] ■ To explain 
this difference, one may recall that in the spin system two 
states (up and down) are degenerate and there is no net 
gradient potential, which drives system from the up state 
to the down state or vice versa. However, in the case of 
DNA, there is a net gradient potential term arising due 
to the attractive interaction between complimentary nu- 
cleotides, which always drives the system from the open 
state to the closed state during each cycle. At high am- 
plitude, we did not get reliable scaling, where all curves 
collapse on a single curve. In single molecule experi- 
ments, measurements have been done at non-equilibrium 
condition. It is always desirable to infer the equilibrium 
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properties of the system from these data to have better 
understanding of the system. In order to do so, gener- 
ally measurements have been taken in the quasi static 
limit [i| so that techniques involved in thermodynam- 
ics can be employed. In recent years, there have been 
considerable work to extract equilibrium properties from 
the non-equilibrium data e.g. Jarzynski equality, which 
relates free energy differences between two equilibrium 
states through non-equilibrium processes [30| . In another 
approach dominant reaction pathway algorithm is devel- 
oped to compute the most probable reaction pathways 
between two equilibrium states [3lJ. Here, we use the 
work theorem to derive the equilibrium path between the 
two states [22]. The underlying assumption behind the 
work theorem relies on the fact that the initial state of 
the system should be in thermal equilibrium. 

Instead of repeating the force cycle C times, we now 
randomly choose any C conformations, which belong to 
equilibrium conformations at that T(F = 0). We fol- 
low the similar protocol described above to reach the fi- 
nal state from the initial state. No attempt is made to 
achieve equilibrium during this process. The total work 
performed on the system going from the zipped state to 
open state (forward path) is w m = — A/J^J^ yi. When 
the applied force decreases (backward path), the work 
done by the system from the open state to the zipped 
state can be written as w\ — A/^J =m j/j. The equi- 
librium distance y^ for the force fk for the forward path 
can be obtained by assigning the weight exp(—/3wk) to 
all forward C paths 22] at that instant fc, which can be 
written as 



J2i ykexp{-/3w k ) 
T,i exp(-Pw k ) 



(3) 
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FIG. 5: Scaling of loop area of hysteresis Ai oop with respect 
to i/°- 5 F ' 5 in the low frequency limit (a) and with respect to 
v~ x F h in the high frequency limit (b). 

Similarly, yu for the reverse path can also be obtained. 

In Fig. 5. we have plotted the simple average of ex- 
tension over many (C = 1000) forward paths as well as 
backward paths (n = 10 4 LD time steps). One can see 
the existence of hysteresis. In this figure, the solid line 
corresponds to the equilibrium path, which is the same, 
whether we start from the zipped or open conformations. 
In our simulation, the data has been sampled over four 
times of the equilibration time. We have used 2 x 10 9 time 
steps out of which the first 5 x 10 s steps have not been 
taken in the averaging. The results are averaged over 
many trajectories, which are almost the same within the 
standard deviation. The weighted average of the exten- 
sion for the forward and backward paths have also been 
depicted in this plot. One can see from these plots that 
the weighted average even for n — 10 4 LD steps is quite 
close to the equilibrium value (5 x 10 s LD steps). We fur- 
ther note that the weighted average of the backward path 
almost overlaps with the equilibrium path. This may be 
because of the fact that in a reverse path, two strands 
of DNA are in the open state and the system can access 
more configurational space. This gives the higher prob- 
ability of choosing the rare conformations, which have 
dominant contributions in Eq. 3. 

Force induced transition in bio-molecules is the sub- 



ject of interest for the last two decades [13j. However, 
there is almost no effort to study the dynamic transition, 
which can be induced by the oscillatory force. For the 
first time, we show the existence of dynamic transition in 
DNA unzipping, which can be observed in vitro as well as 
in other systems e.g. protein unfolding and RNA unfold- 
ing. In low frequency limit, the value of exponents a and 
f3 associated with this transition are in agreement with 
the values obtained for the isotropic spin case. Because 
of the asymmetry in the present case, in high frequency 
limit, the scaling obtained here differs significantly from 
the isotropic spin case. This is because the open state 
is always under the drive towards the zipped state under 
the potential gradient term, which is absent in the spin 
system. We also showed that using the work theorem, 
it is possible to extract the equilibrium properties of the 
system from the non-equilibrium data, which have poten- 
tial application in single molecule force measurements. 
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FIG. 6: The variation of average extension with cyclic force of 
amplitude 0.4. Dotted and dashed lines show the simple aver- 
aging, where as up and down triangles correspond to weighted 
average, of backward and forward paths, respectively. The 
solid line represents the equilibrium / — x curve 
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